close all;
clear all;

Cellsize = 5;       
Nr=1; %number of the runs\trajectories to read

%% Calculation of the patch geometric center evolution in time; to be stored in mcoord = array(number of runs,number of frames, 2); 2=>[xind, yind] indexes of the patch geom center.
for i=1:Nr %in case we are averaging over many runs. 
    %fname=fullfile('.',num2str(i),'/data_Cdc42T_time_course') %for many numbered trajectories
    fname=fullfile('.','/data_Cdc42T_time_course') %for 1 trajectory in the current folder.
    clear cdct;
    cdct=importdata(fname,' ');
    if i==1 %detect sizing and create empty matrices in the begining. Assumption: all runs have the same grid size N and number of frames Nf.
        N=size(cdct,2); % N, size of the grid
        Nf=size(cdct,1)/N; %number of frames scales with N, size of the grid.
        mcoord=zeros(Nr,Nf,2); %mcoord = array(number of runs,number of frames, 2); 2=>[xind, yind] indexes of the patch geom center.
        ddx  = Cellsize*sqrt(pi)/N; %grid scale (micron)
    end
    if mod(size(cdct,1),size(cdct,2))~=0
        display('Some frames were not fully written! The further progression of the script will be erroneous\n');
    else
        for j=0:(Nf-1)
            cdci=cdct(1+j*N:N+j*N,:);
            
            % "naive" scheme
%             [xind,yind]=find(cdci>0.7*max(max(cdci))); 
%             xmean=mean(xind);
%             ymean=mean(yind);

            %This scheme follows a patch geom. center, even if it goes over
            %the periodic boundary. The sudden jump over the box\grid
            %length will be taken into account in the msd routine    
            [xind,yind]=find(cdci>0.7*max(max(cdci))); %the value of the threshold can be varied
            if((max(xind)- min(xind))>N/2) %space between found patches is larger than half of the box->the patch is going over the boundary
                s1=xind(find(xind<N/2)); % bins found on the left side
                s2=xind(find(xind>=N/2)); % bins found on the right side
                if (size(s2)>size(s1)) xmean=mean(vertcat(s1+N,s2)); if (xmean>N) xmean=xmean-N; end %the larger part of the patch is on the left.
                else xmean=mean(vertcat(s1,s2-N)); if (xmean<0) xmean=xmean+N; end %on the right. Note: we have to pick one side anyway, including the case s1=s2.
                end
            else xmean=mean(xind); %the patch is grouped in one part of the box, use the geometric average.
            end
            if((max(yind)- min(yind))>N/2) %space between found patches is larger than half of the box->the patch is going over the boundary
                s1=yind(find(yind<N/2)); % bins found on the top
                s2=yind(find(yind>=N/2)); % bins found on the bottom
                if (size(s2)>size(s1)) ymean=mean(vertcat(s1+N,s2)); if (ymean>N) ymean=ymean-N; end%the larger part of the patch is on the left.
                else ymean=mean(vertcat(s1,s2-N)); if (ymean<0) ymean=ymean+N; end%on the right. Note: we have to pick one side anyway, including the case s1=s2.
                end
            else ymean=mean(yind); %the patch is grouped in one part of the box, use the geometric average.
            end
            mcoord(i,j+1,:)=[xmean, ymean];
        end
    end
end

figure(3);
hold on;
CM = jet(Nr);
for i=1:Nr
    scatter(mcoord(i,:,1),mcoord(i,:,2), 25,CM(i,:));
    plot(mcoord(i,:,1),mcoord(i,:,2),'color',CM(i,:));
    pause;
    xlim([0 N]);
    ylim([0 N]);
    title('Movement of the geometric center of Cdc42T patch','fontsize', 20);
end
save mcoord_1run_norec_shift_1.mat mcoord;

%% M.S.D. calculation: take into account periodic boundary conditions(pbc).
dt=10/60; %delta t in min from the RUN_.....
msd=[];t=[];mstd=[];
for j=1:Nf
    sd_j=[];
    for k=1:Nr
        i=1;
        while (i+j)<=Nf
            
            dx=mcoord(k,i+j,1)-mcoord(k,i,1);
            if (dx<0 && abs(dx)>N/2) %jump over a right boundary to the left side. Minimal image is on the left
                dx=dx+N;
            elseif (dx>N/2) %jump over the left boundary to the right side. Minimal image is on the right.
                dx=dx-N; %sign does not matter, we need only the squared value of the minimal distance in pbc.
            end
            dy=mcoord(k,i+j,2)-mcoord(k,i,2);
            if (dy<0 && abs(dy)>N/2) %jump over a bottom boundary to the top. Minimal image is in the bottom
                dy=dy+N;
            elseif (dy>N/2) %jump over the top boundary to the bottom. Minimal image is on the top.
                dy=dy-N; %sign does not matter, we need only the squared value of the minimal distance in pbc.
            end
            sd_j=[sd_j (dx^2 + dy^2)*ddx^2];
            i=i+1;
        end
    end
    %hist(sd_j,0:0.01:max(sd_j));
    %size(sd_j)
    msd=[msd mean(sd_j)];
    mstd=[mstd std(sd_j)];
    t=[t dt*j]; %time in minutes
end
figure(110);
hold on;
plot(t,msd,'-b', 'LineWidth',3);
%errorbar(t,msd,mstd);
xlim([0 20]);
%ylim([0 6]);
xlabel('Time, (min)','fontsize',20);
ylabel('$$\langle \Delta X^2 \rangle$$, (micron)','interpreter','latex','fontsize',20);
title(['MSD for Cdc42T geom center for ' num2str(t(end)/60,3) ' hr simulation'],'fontsize',22);
save msd_1run_norec_shift_1.mat msd;       